function [states,chains] = sim_markov_chain_fast(P,T,N,s0,V)
                                   

%  PURPOSE:
%  Generates N simulations of length T from a Markov chain with transition 
%  probability matrix P and states V, given an innitial state s0 for each 
%  of the N simulations
%
%  OUTPUTS:
%    states : the realized states for each of N chains over T periods
%    chains : the realized values for each of N chains over T periods
%
%  INPUTS:
%    P  : transition matrix
%    T  : number of periods to simulate
%    N  : number of MCs to simulate
%    s0 : initial states, dimension [N,1]
%    V  : quantity corresponding to each state


% auxiliary objects
cdf = [zeros(1,size(P,1)) ; cumsum(P,2)']; % cumulative distribution for each current state [S,S+1]

% initialize output                      
states      = NaN(T,N); 
states(1,:) = s0; 
  
% simulate
if N > 1
   for t = 1:(T-1) % loop over periods
       half_U        = rand(1,N/2);
       U             = [half_U,1-half_U]; % stratify
       states(t+1,:) = sum(bsxfun(@gt,U,cdf(:,states(t,:))));         
   end
elseif N == 1 
   U = rand(T-1,1);
   for t = 1:(T-1)  % loop over periods
       states(t+1) = sum(U(t)>cdf(:,states(t)));
   end  
end
   
if nargout>1
   chains = V(states);
end








% % auxiliary objects
%   cdf = [zeros(size(P,1),1) , cumsum(P,2)]; % cumulative distribution for each current state [S,S+1]
% 
% % initialize output                      
%   states      = NaN(N,T); 
%   states(:,1) = s0; 
%   
% % simulate
%   if N > 1
%      for t = 1:(T-1); % loop over periods
%          U             = rand(1,N/2);
%          U             = [U,1-U]'; % stratify
%          states(:,t+1) = sum(bsxfun(@gt,U,cdf(states(:,t),:)),2);
%      end
%   elseif N == 1 
%      U = rand(T-1,1);
%      for t = 1:(T-1); % loop over periods
%          states(:,t+1) = sum(U(t)>cdf(states(t),:));
%      end  
%   end
%  
%   
% if nargout>1
%    chains = V(states);
%    if N==1, chains = chains'; states = states'; end;
% end